#### Replication for Tables and Figures in Appendix ####
# Dylan Potts, EUI

#Packages
library(pacman)
p_load(tidyverse, estimatr, stargazer, fixest, fwildclusterboot, jtools, dotwhisker, marginaleffects,
       sjPlot, sjmisc, lfe, modelsummary, car, KableExtra, broom.mixed)

# Load Datasets - these can also be found on dataverse
load("Data/datasets/msr_army_trim.RData")
# Case Study datasets
# first the general dataset
load("Data/datasets/monthly_case.RData")
# Next the monthly_dataset
load("Data/datasets/types_case.RData")
# Donations dataset
load("Data/datasets/donations.RData")

# Table 1

m1 <- felm(desert~  irish + german + english| 0 | 0 | comp_id , data=msr_army_trim)
m2 <- felm(desert~ comp_id + irish + german + english | comp_id  | 0 |comp_id, data=msr_army_trim)
m3 <- felm(desert~ comp_id + enlisted_year + irish + german + english | comp_id + enlisted_year | 0 |comp_id, data=msr_army_trim)
m4 <- felm(desert~ irish*fourties + german*fourties + english*fourties | 0 | 0 |comp_id, data=msr_army_trim)
m5 <- felm(desert~ comp_id + irish*fourties + german*fourties + english*fourties | comp_id  | 0 |comp_id, data=msr_army_trim)
m6 <- felm(desert~ irish*fourties + german*fourties + english*fourties | comp_id + enlisted_year | 0 |comp_id, data=msr_army_trim)


stargazer::stargazer(m1, m2, m3, m4, m5, m6, 
                     keep = c("irish", "german", "english", "irish:fourties",
                              "german:fourties", "english:fourties", "fourties"), 
                     title="Effect of Cohorts on Desertion",
                     style="default",
                     keep.stat=c("n"),
                     add.lines=list(c("Company FEs", "", "Y", "Y", "", "Y","Y" ), c("Enlistment Year FEs", "", "", "Y", "", "","Y")),
                     covariate.labels = c("Irish", "1840s",
                                          "German", "English", "Irish x 1840s", "German x 1840s", "English x 1840s"),
                     order=c(1,3,4,2,5,6),
                     dep.var.caption  = "Probability of Desertion", 
                     dep.var.labels   = "", 
                     notes = "Standard errors clustered by company",
                     header=FALSE, 
                     type="latex", 
                     label = "tab:cohortlin"
)

# Table 2

m1a <- glm(desert~  irish + german + english, data=msr_army_trim, family = 'binomial')
m2a <- glm(desert~ comp_id + irish + german + english, data=msr_army_trim, family = 'binomial')
m3a <- glm(desert~ comp_id + enlisted_year + irish + german + english, data=msr_army_trim, family = 'binomial')
m4a <- glm(desert~ irish*fourties + german*fourties + english*fourties, data=msr_army_trim, family = 'binomial')
m5a <- glm(desert~ comp_id + irish*fourties + german*fourties + english*fourties, data=msr_army_trim, family = 'binomial')
m6a <- glm(desert~ comp_id + enlisted_year + irish*fourties + german*fourties + english*fourties, data=msr_army_trim, family = 'binomial')

clust <- function(m, cluster="comp_id") { return(coef(summary(m, cluster = c(cluster)))[, 2]) }

stargazer::stargazer(m1a, m2a, m3a, m4a, m5a, m6a, 
                     
                     keep = c("irish", "german", "english", "fourties", "irish:fourties",
                              "fourties:german", "fourties:english"), 
                     
                     title="Effect of Cohorts on Desertion (Logistic Regression)",
                     style="default",
                     keep.stat=c("n"),
                     add.lines=list(c("Company FEs", "", "Y", "Y", "", "Y","Y" ), c("Enlistment Year FEs", "", "", "Y", "", "","Y")),
                     covariate.labels = c("Irish", "1840s ", "German","English", "Irish x 1840s", "German x 1840s", "English x 1840s"),
                     dep.var.caption  = "Probability of Desertion", 
                     dep.var.labels   = "", 
                     header=FALSE, 
                     type="latex", 
                     notes.label = "Standard errors clustered by company",
                     label = "tab:cohortlog"
)

# Figure 1

msr_army_trim <- msr_army_trim %>% 
  mutate(thse = ifelse(ybirth>=1837, 1, 0),
         thei = ifelse(ybirth>=1838, 1, 0), 
         thni = ifelse(ybirth>=1839, 1, 0), 
         foon = ifelse(ybirth>=1841, 1, 0), 
         fotw = ifelse(ybirth>=1842, 1, 0), 
         foth = ifelse(ybirth>=1843, 1, 0), 
         fofo = ifelse(ybirth>=1844, 1, 0), 
         thtw = ifelse(ybirth>=1832, 1, 0), 
         twfi = ifelse(ybirth>=1825, 1, 0))

m37 <- feols(desert~ irish*thse | enlisted_year + comp_id, cluster = 'comp_id', data=msr_army_trim)
m38 <- feols(desert~ irish*thei | enlisted_year + comp_id, cluster = 'comp_id', data=msr_army_trim)
m39 <- feols(desert~ irish*thni | enlisted_year + comp_id, cluster = 'comp_id', data=msr_army_trim)
m40 <- feols(desert~ irish*fourties | enlisted_year + comp_id, cluster = 'comp_id', data=msr_army_trim)
m41 <- feols(desert~ irish*foon | enlisted_year + comp_id, cluster = 'comp_id', data=msr_army_trim)
m42 <- feols(desert~ irish*fotw | enlisted_year + comp_id, cluster = 'comp_id', data=msr_army_trim)
m43 <- feols(desert~ irish*foth | enlisted_year + comp_id, cluster = 'comp_id', data=msr_army_trim)
m44 <- feols(desert~ irish*fofo | enlisted_year + comp_id, cluster = 'comp_id', data=msr_army_trim)
m32 <- feols(desert~ irish*thtw | enlisted_year + comp_id, cluster = 'comp_id', data=msr_army_trim)
m25 <- feols(desert~ irish*twfi | enlisted_year + comp_id, cluster = 'comp_id', data=msr_army_trim)


cutoff_coef <- plot_coefs(m37, m38, m39, m40, m41, m42, m43,
                          coefs = c("1837"= "irish:thse", 
                                    "1838" = "irish:thei",
                                    "1839" = "irish:thni", 
                                    "1840" = "irish:fourties", 
                                    "1841" = "irish:foon",
                                    "1842" = "irish:fotw",
                                    "1843" = "irish:foth", 
                                    "1844" = "irish:fofo"),
                          colors = c( "Aquamarine3", "Aquamarine3", "Aquamarine3", "Aquamarine3", "Aquamarine3", "Aquamarine3", "Aquamarine3"),
                          scale = TRUE, 
                          point.shape=FALSE, 
                          point.size=4, 
                          secret_weapon=TRUE
                          #groups = list(pane_1 = c("younger"),pane_2 = c("youngerge"), pane_3 = c("youngeram")))
)

cutoff_coef <- cutoff_coef + theme(legend.position="none") +coord_flip()
cutoff_coef <- cutoff_coef +
  ggtitle("")+
  xlab("Desertion")+
  ylab("Famine Cohort Cutoff") +
  theme(axis.text.x = element_text(size = 17))



# Table 4
famir3l <- glm(desert ~ younger  + famid + enlisted_year, data= (msr_army_trim %>% filter(irish==1) %>% group_by(famid) %>% filter(n()>1) %>%
                                                                   filter(n_distinct(younger) ==2) %>% ungroup()), family = "binomial")

famge3l <- glm(desert ~ youngerge + famid + enlisted_year, data= (msr_army_trim %>% filter(german==1) %>% group_by(famid) %>% filter(n()>1) %>%
                                                                    filter(n_distinct(younger) ==2) %>% ungroup()), family = "binomial")

famam3l <- glm(desert ~ youngeram + famid + enlisted_year, data= (msr_army_trim %>% filter(american==1) %>% group_by(famid) %>% filter(n()>1) %>%
                                                                    filter(n_distinct(younger) ==2) %>% ungroup()), family = "binomial")
                                                                    
                                                             
                                                                    stargazer::stargazer(famir3l, famge3l, famam3l, 
                                                                                         keep = c("younger", "youngerge", "youngeram"), 
                                                                                         type="latex",
                                                                                         title="Effect of Sibling on Birth Order (Logistic)",
                                                                                         style="default",
                                                                                         keep.stat=c("n"),
                                                                                         add.lines=list(c("Family FEs", "Y", "Y", "Y"), c("Enlistment Year FEs", "Y", "Y", "Y")),
                                                                                         covariate.labels = c("Younger Irish", "Younger German", "Younger American"),
                                                                                         dep.var.caption  = "Probability of Desertion", 
                                                                                         dep.var.labels   = "", 
                                                                                         column.labels   = c("Irish Siblings", "German Siblings", "American Siblings"),
                                                                                         column.separate = c(1, 1),
                                                                                         header=FALSE, 
                                                                                         
                                                                                         label="tab:logitsib" )                                                                   


# Figure 2

# irish
famirap1 <- feols(desert ~ younger  | famid + enlisted_year, data= (msr_army_trim %>% filter(irish==1) %>% 
                                                                      add_count(last, name="GSIZE") %>% filter(GSIZE==2) %>% 
                                                                      group_by(famid) %>% filter(n()>1) %>%
                                                                      filter(n_distinct(younger) ==2) %>% ungroup()))
famirap2 <- feols(desert ~ younger  | famid + enlisted_year, data= (msr_army_trim %>% filter(irish==1) %>% 
                                                                      add_count(last, name="GSIZE") %>% filter(GSIZE==2) %>% 
                                                                      group_by(famid) %>% filter(n()>1) %>%
                                                                      filter(n_distinct(younger) ==2) %>% ungroup()))
famirap3 <- feols(desert ~ younger  | famid + enlisted_year, data= (msr_army_trim %>% filter(irish==1) %>% 
                                                                      add_count(last, name="GSIZE") %>% filter(GSIZE==2) %>% 
                                                                      group_by(famid) %>% filter(n()>1) %>%
                                                                      filter(n_distinct(younger) ==2) %>% ungroup()))
# german

famgeap1 <- feols(desert ~ younger | famid + enlisted_year, data= (msr_army_trim %>% filter(german==1) %>%  
                                                                     add_count(last, name="GSIZE") %>% filter(GSIZE==2) %>%
                                                                     group_by(famid) %>% filter(n()>1) %>%
                                                                     filter(n_distinct(younger) ==2) %>% ungroup()))
famgeap2 <- feols(desert ~ younger | famid + enlisted_year, data= (msr_army_trim %>% filter(german==1) %>%  
                                                                     add_count(last, name="GSIZE") %>% filter(GSIZE==2) %>%
                                                                     group_by(famid) %>% filter(n()>1) %>%
                                                                     filter(n_distinct(younger) ==2) %>% ungroup()))
famgeap3 <- feols(desert ~ younger | famid + enlisted_year, data= (msr_army_trim %>% filter(german==1) %>%  
                                                                     add_count(last, name="GSIZE") %>% filter(GSIZE==2) %>%
                                                                     group_by(famid) %>% filter(n()>1) %>%
                                                                     filter(n_distinct(younger) ==2) %>% ungroup()))
# american

famamap1 <- feols(desert ~ younger | famid + enlisted_year, data= (msr_army_trim %>% filter(american==1) %>% 
                                                                     add_count(last, name="GSIZE") %>% filter(GSIZE==2) %>% 
                                                                     group_by(famid) %>% filter(n()>1) %>%
                                                                     filter(n_distinct(younger) == 2) %>% ungroup()))
famamap2 <- feols(desert ~ younger | famid + enlisted_year, data= (msr_army_trim %>% filter(american==1) %>% 
                                                                     add_count(last, name="GSIZE") %>% filter(GSIZE==2) %>% 
                                                                     group_by(famid) %>% filter(n()>1) %>%
                                                                     filter(n_distinct(younger) == 2) %>% ungroup()))
famamap3 <- feols(desert ~ younger | famid + enlisted_year, data= (msr_army_trim %>% filter(american==1) %>% 
                                                                     add_count(last, name="GSIZE") %>% filter(GSIZE==2) %>% 
                                                                     group_by(famid) %>% filter(n()>1) %>% 
                                                                     filter(n_distinct(younger) == 2) %>% ungroup()))



sib_coefap <- plot_coefs(famirap3,famirap2, famirap1,famgeap3,famgeap2,famgeap1,famamap3,famamap2,famamap1,
                         coefs = c("Younger American"= "youngeram", 
                                   "Younger German" = "youngerge",
                                   "Younger Irish" = "younger"),
                         scale = TRUE, 
                         colors = c("darkblue", "blue", "Aquamarine3", "darkblue", "blue", "Aquamarine3", "darkblue", "blue", "Aquamarine3" ), 
                         point.shape=FALSE, 
                         point.size=4
                         #groups = list(pane_1 = c("younger"),pane_2 = c("youngerge"), pane_3 = c("youngeram")))
)

sib_coefap <- sib_coefap + theme(legend.position="none")+coord_flip()
sib_coefap <- sib_coefap +
  ggtitle("Estimated Effect of Being a Younger Sibling on Desertion")+
  xlab("Desertion")+
  ylab("Ethnic Group")     



# Figure 3
m_hei3 <- lm(desert~ comp_id + enlisted_year + rh_cm*fourties, data=(msr_army_trim %>% filter(irish==1)%>% filter(rh_cm > 140)))

irint<- plot_cap(m_hei3, condition = c("rh_cm", "fourties"), vcov = ~comp_id) # base plot
irint <-irint + scale_x_continuous(limits = c(145, 185))+
  ggtitle("")+
  ylab("Probability of Desertion")+
  xlab("Soldier Height (cm)")+
  theme_classic()


# Figure 4
m_hei7 <- lm(desert~ comp_id + enlisted_year + rh_cm*fourties, data=(msr_army_trim %>% filter(english==1)))
engint<- plot_cap(m_hei7, condition = c("rh_cm", "fourties"), vcov = ~comp_id)
engint <- engint + scale_x_continuous(limits = c(148, 182))+
  ggtitle("English")+
  ylab("Probability of Desertion")+
  xlab("Soldier Height")+
  theme_classic()

# Figure 5
m_hei8 <- lm(desert~ comp_id + enlisted_year + rh_cm*fourties, data=(msr_army_trim %>% filter(german==1)))
gerint <- plot_cap(m_hei8, condition = c("rh_cm", "fourties"), vcov = ~comp_id)
gerint <- gerint + scale_x_continuous(limits = c(148, 182))+
  ggtitle("Germans")+
  ylab("Probability of Desertion")+
  xlab("Soldier Height")+
  theme_classic()

# Figure 6
m_hei9 <- lm(desert~ comp_id + enlisted_year + rh_cm*fourties, data=(msr_army_trim %>% filter(american==1)))
amint <- plot_cap(m_hei9, condition = c("rh_cm", "fourties"), vcov = ~comp_id)
amint <-amint + scale_x_continuous(limits = c(148, 182)) +
  ggtitle("Americans")+
  ylab("Probability of Desertion")+
  xlab("Soldier Height")+
  theme_classic()

# Table 5
hma1 <- felm(rh_cm ~ fourties | 0 | 0 | comp_id, data = (msr_army_trim %>%  filter(american ==1)))
hma2 <- felm(rh_cm ~ fourties | comp_id | 0 | comp_id, data = (msr_army_trim %>%  filter(american ==1)))
hma3 <- felm(rh_cm ~ fourties | comp_id + enlisted_year | 0 | comp_id, data = (msr_army_trim %>%  filter(american ==1)))

stargazer::stargazer(hma1, hma2, hma3, 
                     keep = c("fourties"), 
                     type="latex",
                     title="Height Predicted by Age Cohort of American Soldiers",
                     style="default",
                     keep.stat=c("n"),
                     add.lines=list(c("Company FEs", "", "Y","Y" ), c("Enlistment Year FEs","", "","Y")),
                     covariate.labels = c("1840s"),
                     dep.var.caption  = "Height at Enlistment", 
                     dep.var.labels   = "", 
                     notes.label = "Standard errors clustered by company",
                     column.separate = c(1,1),
                     header=FALSE, 
                     font.size = "small",
                     label = "tab:htesta")

# Table 6
m_hei1a <- glm(desert~ rh_cm *fourties, data=(msr_army_trim %>% filter(irish==1) %>% filter(rh_cm > 120)), family='binomial')
m_hei2a <- glm(desert~ comp_id + rh_cm*fourties, data=(msr_army_trim %>% filter(irish==1)%>% filter(rh_cm > 120)), family='binomial')
m_hei3a <- glm(desert~ comp_id + enlisted_year + rh_cm*fourties, data=(msr_army_trim %>% filter(irish==1)%>% filter(rh_cm > 120)), family='binomial')
m_hei4a <- glm(desert~ rh_cm*fourties, data=(msr_army_trim %>% filter(irish==0)%>% filter(rh_cm > 120)), family='binomial')
m_hei5a <- glm(desert~ comp_id + rh_cm*fourties, data=(msr_army_trim %>% filter(irish==0)%>% filter(rh_cm > 120)), family='binomial')
m_hei6a <- glm(desert~ comp_id + enlisted_year + rh_cm*fourties, data=(msr_army_trim %>% filter(irish==0) %>% filter(rh_cm > 120)), family='binomial')

heighttaba <- stargazer(m_hei1a, m_hei2a, m_hei3a, m_hei4a, m_hei5a, m_hei6a, 
                        keep = c("rh_cm:fourties", "rh_cm", "fourties"), 
                        type='latex',
                        title="Soldier Height as a Natural Experiment (Logistic Regression)",
                        style="default",
                        keep.stat=c("n", "AIC"),
                        add.lines=list(c("Company FEs", "", "Y", "Y", "", "Y", "Y"), c("Enlistment Year FEs", "", "", "Y", "", "", "Y")),
                        covariate.labels = c("Soldier Height", "Born in the 1840s", "Soldier Height x Born in the 1840s"),
                        dep.var.caption= "Probability of Desertion", 
                        dep.var.labels= "", 
                        notes.label = "Standard errors clustered by company",
                        column.labels   = c("Irish Soldiers", "All Other Soldiers Placebo"),
                        column.separate = c(3, 3),
                        header=FALSE 
)

# Figure 7

corrwl <- ggplot(data=donat, aes(x=works, y=exposure)) +geom_point() +geom_smooth(method=lm)+
  labs(title="Relationship Between Famine Public Works and Population Change") +
  labs(x="Public Works During Famine", y="Ppulation Change")+ theme_bw()


# Table 7

m_ind60 <- felm(desert~ findex + ybirth + manual| enlisted_year + comp_id + province| 0 | comp_id, data=(msr_army_trim %>% filter(share > 0.6)))
m_ind50 <- felm(desert~ findex + ybirth + manual| enlisted_year + comp_id + province| 0 | comp_id, data=msr_army_trim)
m_ind70 <- felm(desert~ findex + ybirth + manual| enlisted_year + comp_id + province| 0 | comp_id, data=(msr_army_trim %>% filter(share > 0.7)))
m_ind80 <- felm(desert~ findex + ybirth + manual| enlisted_year + comp_id + province| 0 | comp_id, data=(msr_army_trim %>% filter(share > 0.8)))
m_ind90 <- felm(desert~ findex + ybirth + manual| enlisted_year + comp_id + province| 0 | comp_id, data=(msr_army_trim %>% filter(share > 0.9)))

stargazer::stargazer(m_ind50, m_ind60, m_ind70, m_ind80, m_ind90, 
                     keep = c("findex"), 
                     type="latex",
                     title="Varying Fixed Effect Construction",
                     style="default",
                     keep.stat=c("n"),
                     add.lines=list(c("Province Confidence percent","50","60","70","80","90"),c("Company FEs","Y","Y","Y","Y","Y"), c("Enlistment Year FEs","Y","Y","Y","Y","Y")),
                     covariate.labels = c("Famine Index"),
                     dep.var.caption  = "Probability of Desertion", 
                     dep.var.labels = "",
                     notes.label = "Standard errors clustered by company",
                     header=FALSE, 
                     font.size = "small",
                     column.separate = c(1, 1),
                     label = "tab:originsp")


# Figure 8
donations <- ggplot(data=donat, aes(x=exposure, y=Donation)) +geom_point() +geom_smooth(method=lm)+
  labs(title="Relationship Between Famine Exposure and Donations", subtitle = "Irish Soldiers' Charity to Food Relief in Ireland (1863)") +
  labs(x="Population Change")+ theme_bw()


# Figure 9
works <- ggplot(data=donat, aes(x=log(1+works), y=Donation)) +geom_point() +geom_smooth(method=lm)+
  labs(title="Relationship Between Famine Public Works and Donations", subtitle = "Irish Soldiers' Charity to Food Relief in Ireland (1863)") +
  labs(x="Public Works During Famine")+ theme_bw()

# Figure 10
losshei <- ggplot(data=(msr_army_trim), aes(x=exposure, y=rh_cm)) +geom_point() +geom_smooth(method=lm)+
  labs(title="Relationship Between Famine Population Change and Adult Height") +
  labs(x="Population Change", y="Adult Height")+ theme_bw()

# Figure 11
workshei <- ggplot(data=(msr_army_trim  %>% filter (irish==1) %>% filter(rh_cm >130)), aes(x=works, y=rh_cm)) +geom_point() +geom_smooth(method=lm)+
  labs(title="Relationship Between Famine Public Works and Adult Height") +
  labs(x="Public Works During Famine", y="Adult Height")+ theme_bw()

# Table 8
msr_army_trim <- msr_army_trim %>% 
  mutate(exposureneg = (-1*exposure)) %>% 
  mutate(works = log(1+works))

m_pop1h <- felm(rh_cm~ exposureneg + manual + ybirth | 0 | 0 | comp_id, data=(msr_army_trim %>%  filter(irish==1)))
m_wor1h <- felm(rh_cm~ log(1+works) + manual + ybirth| 0 | 0 | comp_id, data=(msr_army_trim %>%  filter(irish==1)))
m_fin1h <- felm(rh_cm~ findex + manual + ybirth| 0 | 0 | comp_id, data=(msr_army_trim %>%  filter(irish==1)))

stargazer::stargazer(m_pop1h, m_wor1h, m_fin1h,
                     keep = c("exposureneg", "works", "findex"), 
                     type="latex",
                     title="Effect of Famine Death Rates and Public Works on Adult Heights",
                     style="default",
                     keep.stat=c("n"),
                     covariate.labels = c("Population Loss", "Log Public Works", "Famine Index"),
                     dep.var.caption  = "Adult Height", 
                     dep.var.labels   = "", 
                     notes = c("Standard errors clustered by company", 
                               "Control for soldier birth year and Manual Occupation"),
                     header=FALSE, 
                     label = "tab:originx" 
                     
)
# Table 9
msr_army_trim <- msr_army_trim %>% 
  mutate(exposureneg = (-1*exposure))

m_pop1a <- glm(desert~ exposureneg + manual + ybirth, data=(msr_army_trim %>%  filter(irish==1)), family = 'binomial')
m_pop2a <- glm(desert~ comp_id + exposureneg + manual + ybirth, data=(msr_army_trim %>%  filter(irish==1)), family = 'binomial')
m_pop3a <- glm(desert~ comp_id + enlisted_year + exposureneg + manual + ybirth, data=(msr_army_trim %>%  filter(irish==1)), family = 'binomial')
m_wor1a <- glm(desert~ log(1+works) + manual + ybirth, data=(msr_army_trim %>%  filter(irish==1)), family = 'binomial')
m_wor2a <- glm(desert~ comp_id + log(1+works) + manual + ybirth, data=(msr_army_trim %>%  filter(irish==1)), family = 'binomial')
m_wor3a <- glm(desert~ comp_id + enlisted_year + log(1+works) + manual + ybirth, data=(msr_army_trim %>%  filter(irish==1)), family = 'binomial')


stargazer::stargazer(m_pop1a, m_pop2a, m_pop3a, m_wor1a, m_wor2a, m_wor3a,  
                     
                     keep = c("exposureneg", "works"), 
                     type="latex",
                     title="Effect of Famine Death Rates on Cohorts of Irish Desertion",
                     style="default",
                     keep.stat=c("n"),
                     add.lines=list(c("Company FEs", "", "Y", "Y", "", "Y", "Y"), 
                                    c("Enlistment Year FEs", "", "", "Y", "", "", "Y")),
                     covariate.labels = c("Population Loss", "Log Public Works"),
                     dep.var.caption  = "Probability of Desertion", 
                     dep.var.labels   = "", 
                     notes = c("Standard errors clustered by company", 
                               "Control for soldier birth year and Manual Occupation"),
                     notes.align = "l",
                     header=FALSE, 
                     label = "tab:deathlogit" 
                     
)

# Table 10
mba <- msr_army_trim %>% filter(battles!=0) %>% 
  glm(desert~ comp_id + enlisted_year + irish*fourties + german*fourties + english*fourties, data=.,
      family='binomial')
mnba <- msr_army_trim %>% filter(battles==0) %>%
  glm(desert~ comp_id + enlisted_year + irish*fourties + german*fourties + english*fourties, data=.,
      family = 'binomial')

splittaba <- stargazer::stargazer(mba, mnba, 
                                  keep = c("irish", "fourties", "irish:fourties"), 
                                  type="latex",
                                  title="Effects of Cohorts When Facing Major Battles (Logistic)",
                                  style="default",
                                  keep.stat=c("n", "AIC"),
                                  add.lines=list(c("Company FEs", "Y","Y" ), c("Enlistment Year FEs", "Y","Y")),
                                  covariate.labels = c("Irish", "German ", "English","1840s", "Irish x 1840s", "German x 1840s", "English x 1840s"),
                                  dep.var.caption  = "Probability of Desertion", 
                                  dep.var.labels   = "", 
                                  notes.label = "Standard errors clustered by company",
                                  column.labels   = c( "At Least One Battle", "No Major Battles" ),
                                  column.separate = c(1,1),
                                  header=FALSE)

# Table 11

m_ment1 <-lm(mental ~ irish + german + english, data=msr_army_trim)
m_ment2 <- lm(mental ~ comp_id +irish + german + english, data=msr_army_trim)
m_ment3 <- lm(mental ~ comp_id + enlisted_year + irish + german + english, data=msr_army_trim)


mentaltab <- stargazer::stargazer(m_ment1, m_ment2, m_ment3, 
                                  keep = c("irish", "german", "english"), 
                                  type="latex",
                                  title="Mental Health Problems Among Different Immigrant Groups",
                                  style="default",
                                  keep.stat=c("n"),
                                  add.lines=list(c("Company FEs", "", "Y","Y" ), c("Enlistment Year FEs","", "","Y")),
                                  covariate.labels = c("Irish", "German ", "English"),
                                  dep.var.caption  = "Probability of Desertion", 
                                  dep.var.labels   = "", 
                                  notes.label = "Standard errors clustered by company",
                                  
                                  column.separate = c(1,1),
                                  header=FALSE)
# Figure 13
ir_comp_a <- lm(desert~ enlisted_year + irish*comp_ir, data=msr_army_trim)
ir_comp_moda <- plot_cap(ir_comp_a, condition = c("comp_ir", "irish"))
ir_comp_moda <- ir_comp_mod + xlab("Share of a Company Which is Irish") +
  ylab("Predicted Desertion Rate") +
  theme_classic()

# Table 12

# Create the values for average arrival time by cohort and nationalist
msr_army_trim$rentry_yr <- as.numeric(msr_army_trim$rentry_yr)

msr_ir_pre40 <- msr_army_trim %>% filter(irish==1, fourties ==0)
irish_pre40s <- mean(msr_ir_pre40$rentry_yr, na.rm = TRUE)

msr_ir_post40 <- msr_army_trim %>% filter(irish==1, fourties ==1)
irish_post40s <- mean(msr_ir_post40$rentry_yr, na.rm = TRUE)

msr_ge_pre40 <- msr_army_trim %>% filter(german==1, fourties ==0)
ger_pre40s <- mean(msr_ge_pre40$rentry_yr, na.rm = TRUE)

msr_ge_post40 <- msr_army_trim %>% filter(german==1, fourties ==1)
ger_post40s <- mean(msr_ge_post40$rentry_yr, na.rm = TRUE)

msr_eng_pre40 <- msr_army_trim %>% filter(english==1, fourties ==0)
eng_pre40s <- mean(msr_eng_pre40$rentry_yr, na.rm = TRUE)

msr_eng_post40 <- msr_army_trim %>% filter(english==1, fourties ==1)
eng_post40s <- mean(msr_eng_post40$rentry_yr, na.rm = TRUE)


# Table 13
msr_army_trim <- msr_army_trim %>% ungroup() %>% 
  mutate(arrival = rentry_yr)
msr_army_trim$arrival <- as.numeric(msr_army_trim$arrival)

m_entry1 <- lm(desert~ arrival, data=msr_army_trim)
m_entry2 <- lm(desert~ arrival, data=(msr_army_trim %>% filter(irish==1)))

stargazer::stargazer(m_entry1, m_entry2,
                     type="latex",
                     title="Effect of Immigrant Arrival Time on Desertion",
                     style="default",
                     keep.stat=c("n"),
                     covariate.labels = c("Arrival Year"),
                     column.labels   = c( "All Immigrants", "Irish Only" ),
                     column.separate = c(1,1),
                     dep.var.caption  = "Probability of Desertion", 
                     dep.var.labels   = "", 
                     header=FALSE, 
                     notes.label = "Standard errors clustered by company",
                     label = "tab:entry" 
                     
)

# Table 14
msr_army_trim <- msr_army_trim %>% 
  mutate(wound=ifelse(milcact1=="W", 1, 0))

m_weak1 <- lm(wound~ rh_cm * fourties, data=(msr_army_trim %>% filter(irish==1)))
m_weak2 <- lm(wound~ comp_id + rh_cm * fourties, data=(msr_army_trim %>% filter(irish==1)))
m_weak3 <- lm(wound~ comp_id + enlisted_year + rh_cm * fourties, data=(msr_army_trim %>% filter(irish==1)))
m_weak4 <- lm(wound~ rh_cm * fourties, data=(msr_army_trim %>% filter(irish==0)))
m_weak5 <- lm(wound~ comp_id + rh_cm * fourties, data=(msr_army_trim %>% filter(irish==0)))
m_weak6 <- lm(wound~ comp_id + enlisted_year + rh_cm * fourties, data=(msr_army_trim %>% filter(irish==0)))



stargazer::stargazer(m_weak1, m_weak2, m_weak3,m_weak4, m_weak5, m_weak6, 
                     se=list(clust(m_weak1), clust(m_weak2), clust(m_weak3), clust(m_weak4), clust(m_weak5), clust(m_weak6)),
                     keep = c("rh_cm:fourties", "rh_cm", "fourties"), 
                     type="latex",
                     title="Soldier Height as a Predictor of Wounding (Physical Disadvantage)",
                     style="default",
                     keep.stat=c("n"),
                     add.lines=list(c("Company FEs", "", "Y", "Y", "", "Y", "Y"), c("Enlistment Year FEs", "", "", "Y", "", "", "Y")),
                     covariate.labels = c("Height", "Born in the 1840s", "Height x Born in the 1840s"),
                     dep.var.caption  = "Probability of Wounding", 
                     dep.var.labels   = "", 
                     column.labels   = c("Irish Soldiers", "All Other Soldiers Placebo"),
                     column.separate = c(3, 3),
                     header=FALSE, 
                     notes = "Standard Errors Clustered at the Company Level", 
                     label="tab:weakt"
)

# Table 15
m_role1 <- lm(battles~ rh_cm * fourties, data=(msr_army_trim %>% filter(irish==1)))
m_role2 <- lm(battles~ comp_id + rh_cm * fourties, data=(msr_army_trim %>% filter(irish==1)))
m_role3 <- lm(battles~ comp_id + enlisted_year + rh_cm * fourties, data=(msr_army_trim %>% filter(irish==1)))
m_role4 <- lm(battles~ rh_cm * fourties, data=(msr_army_trim %>% filter(irish==0)))
m_role5 <- lm(battles~ comp_id + rh_cm * fourties, data=(msr_army_trim %>% filter(irish==0)))
m_role6 <- lm(battles~ comp_id + enlisted_year + rh_cm * fourties, data=(msr_army_trim %>% filter(irish==0)))



stargazer::stargazer(m_role1, m_role2, m_role3,m_role4, m_role5, m_role6, 
                     se=list(clust(m_role1), clust(m_role2), clust(m_role3), clust(m_role4), clust(m_role5), clust(m_role6)),
                     keep = c("rh_cm:fourties", "rh_cm", "fourties"), 
                     type="latex",
                     title="Soldier Height as a Predictor of Facing Battle - Differential Roles in Service?",
                     style="default",
                     keep.stat=c("n"),
                     add.lines=list(c("Company FEs", "", "Y", "Y", "", "Y", "Y"), c("Enlistment Year FEs", "", "", "Y", "", "", "Y")),
                     covariate.labels = c("Height", "Born in the 1840s", "Height x Born in the 1840s"),
                     dep.var.caption  = "Probability of Major Battle", 
                     dep.var.labels   = "", 
                     column.labels   = c("Irish Soldiers", "All Other Soldiers Placebo"),
                     column.separate = c(3, 3),
                     header=FALSE, 
                     notes = "Standard Errors Clustered at the Company Level", 
                     label="tab:role"
)

# Table 16
lmc <- felm(cook~ irish*fourties | enlisted_year + comp_id| 0 | comp_id, data=msr_army_trim)
lmm <- felm(medical~ irish*fourties | enlisted_year + comp_id| 0 | comp_id, data=msr_army_trim)
lmp <- felm(picket~ irish*fourties | enlisted_year + comp_id| 0 | comp_id, data=msr_army_trim)
lmt <- felm(teamster~ irish*fourties | enlisted_year + comp_id| 0 | comp_id, data=msr_army_trim)
lmpr <- felm(prov~ irish*fourties | enlisted_year + comp_id| 0 | comp_id , data=msr_army_trim)
lmpi <- felm(pioneer~ irish*fourties | enlisted_year + comp_id| 0 | comp_id, data=msr_army_trim)

lmtd <- felm(desert~ teamster | enlisted_year + comp_id| 0 | comp_id, data=msr_army_trim)


stargazer::stargazer(lmc, lmm, lmp, lmt, lmpr, lmpi, lmtd, 
                     
                     keep = c("irish", "fourties", "irish:fourties", "teamster"), 
                     type="latex",
                     title="Soldier Cohort Predicting Roles in the Army?",
                     style="default",
                     keep.stat=c("n"),
                     add.lines=list(c("Company FEs", "Y", "Y", "Y", "Y", "Y", "Y", "Y"), 
                                    c("Enlistment Year FEs", "Y", "Y", "Y", "Y", "Y", "Y", "Y")),
                     covariate.labels = c("Irish", "Born in the 1840s", "Irish x Born in the 1840s"),
                     dep.var.caption  = "", 
                     dep.var.labels   = "", 
                     column.labels   = c("Cook", "Medical", "Picket", "Teamster", "Provost", "Pioneer", "Desertion"),
                     column.separate = c(1, 1),
                     header=FALSE, 
                     notes = "Standard Errors Clustered at the Company Level", 
                     font.size = "small",
                     label="tab:roles"
)















